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One of the fundamental problems of solar magnetism is understanding the formation 
of long-lived compact magnetic structures, such as sunspots and pores, in the highly 
turbulent upper layer of the solar convective zone. Observations show that the process of 
accumulation of fragmented magnetic fields emerging from the interior into these compact 
structures occurs in the upper layers of the convection zone. In this study we investigate 
the effects of turbulent convection on formation of large-scale inhomogeneous magnetic 
structures by means of Large-Eddy Simulation (LES) for convection in solar-type stars. 
The main idea of this study is the implementation of a new subgrid-scale model for the 
effective Lorentz force in a three-dimensional nonlinear radiative magnetohydrodynamics 
(MHD) code developed for simulating the upper solar convection zone and lower atmo- 
sphere. To this end we derived the energy budget equations, which include the effects 
of the subgrid-scale turbulence on the Lorentz-force, and implemented the new subgrid- 
scale turbulence model (TELF-Model) in a three-dimensional nonlinear MHD LES code. 
Using imposed initial vertical and horizontal uniform magnetic fields in LES with the 
TELF-Model, we have shown that the magnetic flux tubes formation is started when 
the initial mean magnetic field is larger than a threshold value (about 100 G). This is 
in agreement with the theoretical studies by Rogachevskii & Kleeorin (2007). We have 
determined the vertical profiles of the velocity and magnetic fluctuations, total MHD 
energy and anisotropy of turbulent magneto-convection, kinetic and current and cross 
helicities. 



1. Introduction 

Significant progress toward quantitative understanding of the origin and nature of solar 
activity requires large-scale, integrated modeling of the physical conditions in subsurface 
layers of the Sun, where the magnetic field is generated and organized. Realistic MHD 
numerical simulations of subsurface flows and magnetic structures have become achiev- 
able because of the development of fast supercomputer systems and efficient parallel 
computer codes. In conjunction with modern subgrid-scale turbulence models, these nu- 
merical simulations provide an important tool for investigating the complicated physics 
of the upper turbulent convective boundary layer of the Sun. The dynamics of this layer 
is particularly critical for understanding the mechanism of formation of solar magnetic 
structures and predicting solar activity. Sunspots and pores represent one of the oldest 
and most intriguing problem of solar magnetism. Despite the long history of observa- 
tional and theoretical investigations, the mechanism of their formation is still an open 
question. 
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The goal of this study is to investigate the effects of turbulent convection on the for- 
mation of large-scale inhomogeneous magnetic structures by means of Large-Eddy Sim- 
ulation (LES) for convection in solar-type stars. The key idea of this study is the imple- 
mentation of a new subgrid-scale model for effective Lorentz force in a three-dimensional, 
nonlinear radiative MHD code developed for simulating the upper solar convection zone 
and lower atmosphere (Jacoutot et al. 2008a, b; Kitiashvili et al. 2009, 2010). This new 
subgrid-scale model for the effective Lorentz force is based on 

• theoretical findings of the effects of turbulence on the mean Lorentz force (Kleeorin & 
Rogachevskii 1994; Kleeorin et al. 1996; Rogachevskii & Kleeorin 2007), 

• direct numerical simulations (DNS) of this effect for (i) forced stratified turbulence 
(Brandenburg et al. 2010a, b) and (ii) turbulent convection (Kapyla et al. 2010), 

• mean-field numerical modelling (Brandenburg et al. 2010a, b; Kapyla et al. 2010) based 
on parameterizations both of analytic formulae by Rogachevskii & Kleeorin (2007) and 
the results of the DNS. 

The physics of the effects of turbulence on the mean Lorentz force is as follows. The 
combined effects of the turbulent Reynolds and Maxwell stress tensors can lead to a 
local reduction of the total turbulent pressure and hence to the possibility of self-induced 
concentrations of large-scale magnetic fields (Rogachevskii & Kleeorin 2007). Such a 
process may play an important role in the formation of sunspots and active regions in 
the Sun. DNS of forced stratified turbulence and turbulent convection have demonstrated 
that an imposed mean magnetic field results in a decrease of the effective magnetic 
pressure (Brandenburg et al. 2010a, b; Kapyla et al. 2010). This phenomenon is quantified 
by determining the relevant contributions that relate the sum of the turbulent Reynolds 
and Maxwell stresses with the Maxwell stress of the mean magnetic field. 

Using such a parameterization, Brandenburg et al. (2010a, b) and Kapyla et al. (2010) 
have shown by means of the two-dimensional and three-dimensional mean-field numerical 
modelling that an isentropic (or isothermal) density stratified layer becomes unstable in 
the presence of a uniform imposed magnetic field. This large-scale instability results in 
the formation of loop-like magnetic structures that are concentrated at the top of the 
stratified layer. In three dimensions these structures resemble the appearance of bipolar 
magnetic regions in the Sun. The results of DNS and mean-field numerical modelling 
(Brandenburg et al. 2010a,b; Kapyla et al. 2010) are in good agreement with the theo- 
retical predictions by Rogachevskii & Kleeorin (2007). 

In this study we investigate formation of inhomogeneous magnetic structures using a 
three-dimensional radiative MHD code, "SolarBox," developed by A. Wray at NASA Ad- 
vanced Supercomputing Division (Jacoutot et al. 2008a, b; Kitiashvili et al. 2009, 2010). 
This code takes into account several physical phenomena: compressible fluid flow in a 
highly stratified medium; three-dimensional multigroup radiative energy transfer between 
the fluid elements; a real-gas equation of state, ionization and excitation of all abundant 
species; and magnetic effects. The code is built for simulations of fluid flows in top layers 
of the convective zone and the low atmosphere, in the rectangular geometry. An impor- 
tant feature of this code is implementation of various subgrid-scale turbulence models, 
e.g., the most widely used Smagorinsky model (Smagorinsky 1963) in the compressible 
formulation (Moin et al. 1991; Germano et al. 1991). The turbulent electrical conductiv- 
ity in the code is determined by using the extension of the Smagorinsky model to the 
MHD case (Theobald et al. 1994). The code is based on the LES approach, and solves 
the grid-cell-averaged equations of the conservation of mass, momentum, and energy. We 
modify equations of this LES code by implementing a new subgrid-scale model for effec- 
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tive Lorentz force, and study the dynamics of magnetic structure formation from initial 
horizontal and vertical uniform magnetic fields in a turbulent convection. 

Most of numerical simulations on magnetic flux emerging (see, e.g., Rempel et al. 
2009, Stein et al. 2010) have been done using initial conditions with already existing 
strongly inhomogeneous large-scale magnetic field. On the other hand, formations of 
the large-scale inhomogeneous magnetic structures from initial uniform horizontal or 
vertical magnetic fields in a small-scale turbulent convection have been studied in LES 
by Kitiashvili et al. (2010) and in DNS by Brandenburg et al. (2010a), Brandenburg et 
al. (2010b) and Kapyla et al. (2010). 



2. Derivation of budget equation for energy of mean fields in small-scale MHD 
turbulence 

The effect of the MHD turbulence on the mean (large-scale) magnetic force is described 
by the turbulent magnetic coefficients Q p = 1 — q p and Q s = 1 — q s , so that the mean 
Lorentz force reads 



re// 



Q PB 2 ) +(B- V)^B, (2.1) 



87T / 47T 



(see Rogachevskii & Kleeorin 2007), where B is the mean magnetic field, and the functions 
q p and q s determine the effect of small-scale MHD turbulence on the mean Lorentz force. 
Thus, the equations for the large-scale fields have the following form: 

p(%- + u v)u = -v( P +^B 2 ) +(B- V)^B + F v + F ext , (2.2) 



^8t J V »7T J ^ MtT 

— = Vx[uxB- jjVxB], (2.3) 

^ + V-(pu) = 0, (2.4) 

where u is the mean velocity, p = pk+p T , Pk is the mean fluid pressure, p T is the turbulent 
pressure, F ext is the external force (for example, the gravitational force F ext — pg), 
g is the free-fall acceleration, F„ is the dissipation force owing to the molecular and 
turbulent viscosities, and 77 is the magnetic diffusion (molecular plus turbulent). The 
turbulent diamagnetic (or paramagnetic) drift velocity and the a effect are not included 
in Eq. (I2.3P because both a and the turbulent diamagnetic drift velocity are much smaller 
than the Alfven velocity for the range investigated. Next, we use the following identity: 



^l+u-vW^^ + V,^,,,,). 12.5.1 



that allows us to rewrite Eq. (|2.2p in the following form: 
d(pui) 



n/ , V i (p« i « J -) = -V i lp+^B' , J +(B- V)^B i + [F tf -rP est ]i. (2.' 



In derivation of Eq. (|2.5I) we used Eq. (|2.4p . 

Let us derive the budget equation for energy of the mean fields. We multiply Eq. (|2.6p 
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by the velocity u, Eq. (|2.3|l by (Q s /47r)B, and add them. The result is given by 

Q 
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where £>t is the density of the power released by the turbulent viscosity and the turbu- 
lent magnetic diffusion, D m is the density of the power released through the molecular 
dissipation. Equation (|2.7p is the budget equation for the total energy of the large-scale 



flow and magnetic field pu 2 /2 + Q s B 2 /8tt. The term oc (Q p - Q s ) B 2 (V • u) in RHS of 



Eq. (|2.7p is the work of turbulent magnetic stresses and pressure; the last term in RHS 
of this equation, which depends on time and spatial derivatives of Q s , is the work of the 
pondermotor force. For instance, the term oc dQ s /dt is similar to that in the classical 
electrodynamics when the magnetic permeability depends on time. In the equation for 
the total turbulent energy these terms have opposite signs [see Eq. (|2.12p ], so that the 
equation for the total energy, including the energy of the mean fields, the total (kinetic 
and magnetic) turbulent energy, the potential energy and the internal energy, does not 
contains these terms. In derivation of Eq. (|2.7[) we used the following identities: 
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Using Eq. (|2.7p we obtain budget equation for the total energy (the sum of the energy of 
the mean fields, the potential energy and the internal energy). In particular, the equation 
for the entropy reads 



f + (u.V)S| = / 



pW 7 
~tr~ 



(2.10) 



where S = ln(pp _7 )/7 is the entropy, 7 is the ratio of the specific heats, / is the external 
source of the thermal energy, Wt is the density of the total (kinetic -I- magnetic) energy 
of the MHD turbulence, tr is the characteristic time of the dissipation of the turbulent 
energy into thermal energy (that is of the order of characteristic turbulent time), and 
$ is the total thermal flux. An expression for the internal energy e is determined by 
the first principle of thermodynamics de = TdS + (pk/p 2 )dp. Using this equation and 
Eqs. (|2.4I) . (|2.10p we get the following budget equation for the internal energy: 



pWj 



D m -pfeV-u - V-[$ + upe]. 



(2.11) 



The budget equation for the turbulent energy Wt reads: 
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Then the budget equation for the total energy after taking into account the MHD tur- 
bulence has the following form: 

d ( *ou 2 B 2 \ 

di\^~ + Qs 8^ + pe ) = - v -q + ( F -f u ) + / + / r- ( 2 -!3) 

The total energy flux q is given by 

q = pu(^+e+P) +^Bx(uxB) + u 9p ; ft B 2 +$. (2.14) 

\ 2 p J 47T 87T 

If Qp 7^ Q s , the MHD turbulence produces additional work. It is converted into the 
energy of the large-scale flow and magnetic field even in the absence of dissipation. 
Decrease of the elasticity of the large-scale magnetic field caused by the generation of 
the magnetic fluctuations in the MHD turbulence is essential for systems with the large 
magnetic Reynolds numbers. In this case the developed MHD turbulence can give rise 
to negative Q p . Such values of Q p (i.e., negative) result in excitations of the large-scale 
MHD instabilities, which draw free energy from the turbulent motions. 



3. Governing equations in the TELF-Model and the results of numerical 
simulations 

We use the following compressible MHD conservation equations for mass, momentum, 
energy and magnetic flux: 



^ + V l {pu i ) = 0, 



d(p Ui) 



dt 
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nj = -2C S A 2 |5| [Sij - - (V fc u fc ) 5ij\ + - C c A 2 |S| 2 % 



(3.5) 
(3.6) 
(3.7) 



T) is the magnetic diffusivity due to electrical conductivity of the fluid, n T = u & A/3 is the 
turbulent magnetic diffusivity, k is the molecular heat diffusivity, k t is the turbulent heat 

1/2 

diffusivity, u is the dynamic viscosity, w A is the rms velocity, u A = A [S 2 +g \V Z T\/T] , 
at the scale A (which is the computational grid step size), fa = Bi/B, Zi is the vertical 
unit vector directed opposite to the gravity field, and <fi is the angle between the vertical 
unit vector z and the magnetic field B. 
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Figure 1. Instantaneous distribution of a magnetic field, taking into account the effect of 
turbulent convection on large-scale Lorentz force obtained in LES with the TELF-Model for the 
case of vertical initial magnetic field 10 G. 

The formulae for the nonlinear functions, q p (/3), q s (/3), <7c(/3), and g a (/3), are given 
in Rogachevskii & Kleeorin (2007), where f3 = B/B eq is the non-dimensional magnetic 
field, B cq — y/pu A is the equipartition field strength based on rms velocity u A , and 
the magnetic Reynolds number at this scale Rm A = u A A/?/. In LES discussed in the 
present study we used the following fitting formulae for these nonlinear functions based 
on parameterizations both of analytic formulae by Rogachevskii & Kleeorin (2007) and 
the results of the DNS (Brandenburg et al. 2010a,b; Kapyla, et al. 2010): 



where the parameters can be chosen as C\ = 50 and G<x = 25. 

The numerical simulation results are obtained for the computational domain of 6.4 x 
6.4 x 5.5 Mm with the grid sizes 50 x 50 x 43 km (256 3 mesh points). For details see 
the numerical set-up of Kitiashvili et al. (2010). The domain includes a top, 5 Mm - 
deep layer of the convective zone, and, 0.5 Mm - deep layer of the low atmosphere. 
The lateral boundary conditions arc periodic, and the top and bottom boundaries are 
closed to mass, momentum, and energy transfer (apart from radiative losses). There is 
no net magnetic field lost or gained through the boundaries, i.e., the integral over volume 
for each component of the field is constant. The initial uniform vertical and horizontal 
magnetic fields of various magnitudes (1, 10, and 100 G) were superimposed on the fully 
developed granular convection. The computation runs were up to 3 h of solar time. 

The results of the numerical simulations are given below, whereby we show the in- 
stantaneous distributions of magnetic field obtained in LES, with the TELF-model, that 
takes into account the effect of subgrid turbulent convection on large-scale Lorentz force. 
We considered two cases of the initial vertical (see Figs.[T][2]) and horizontal (see Figs. [3]|4]) 
uniform magnetic field. These figures show that the magnetic flux tubes formation in a 
fully developed turbulent convection depends on the magnitude of the initial magnetic 
field (e.g., when the initial uniform magnetic field is about 100 G, the magnetic flux tubes 
are formed; see Figs. [2] and 2]). This implies that the process of the magnetic flux tubes 
formation is a threshold effect. 

In this study we determined the different statistical characteristic of the turbulent 




cq 



(3.8) 
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Figure 2. Instantaneous distribution of a magnetic field obtained in LES with the 
TELF-Model for the case of vertical initial magnetic field 100 G. 




Figure 3. Instantaneous distribution of a magnetic field obtained in LES with the 
TELF-Model for the case of horizontal initial magnetic field 10 G. 




Figure 4. Instantaneous distribution of a magnetic field obtained in LES with the 
TELF-Model for the case of horizontal initial magnetic field 100 G. 
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Figure 5. Vertical distribution of the turbulent magnetohydrodynamic energy 
p(u 2 )/2 + (b 2 )/87r for the case of vertical initial magnetic fields of different magnitudes. 



convection and formed magnetic field, e.g., we find the vertical profiles of the velocity 
and magnetic fluctuations, MHD energy, anisotropy of turbulent magneto-convection, 
and kinetic, current and cross helicities. In particular, in Figs. [5][7]we show the vertical 
distributions of the turbulent magnetohydrodynamic energy p(u 2 ) /2 + (b 2 ) /8ir and the 
anisotropy parameters of velocity fluctuations, (u 2 )/(u 2 or ), and magnetic fluctuations, 
(b 2 )/ (b^ or ) . Owing to a strong fluid density contrast, the main turbulent energy is located 
at the low part of the convective layer. On the other hand, the anisotropy parameters 
of velocity and magnetic fluctuations are distributed nearly uniformly in the convective 
layer. 



4. Conclusions 

In this study we derived the energy budget equations that include the effects of the 
small-scale (e.g., subgrid-scale) turbulence on the Lorentz-force for modelling of formation 
of the magnetic self-organized structures in the realistic solar conditions. We implemented 
the new subgrid-scale turbulence model (TELF-Model) in a three-dimensional nonlinear 
MHD LES code and performed LES simulations with the TELF-Model with imposed 
initial vertical and horizontal uniform magnetic fields. We have shown that the magnetic 
flux tubes formation is driven when the initial mean magnetic field is larger than a 
threshold value, in agreement with the theoretical studies by Rogachevskii & Kleeorin 
(2007). They demonstrated that the magnetic flux tubes formation is caused by the large- 
scale MHD instability and that a threshold in the initial magnetic field is required for the 
excitation of this instability. We also determined the vertical profiles of the velocity and 
magnetic fluctuations, MHD energy, anisotropy of turbulent magneto-convection, kinetic 
and current and cross helicities. 

In the future study of magnetic flux concentrations one need to determine in LES the 
turbulent Reynolds and Maxwell stress tensors in turbulent convection as a function of 
the initial magnetic field imposed on a turbulent convection. From the latter, one will 
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Figure 6. Vertical distribution of the anisotropy parameter of a turbulent velocity field 
(Uz)/{uj; or } for the case of vertical initial magnetic fields of different magnitudes. 
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Figure 7. Vertical distribution of the anisotropy parameter of a turbulent magnetic field 
(h1)/(b'h OI ) for the case of vertical initial magnetic fields of different magnitudes. 



be able to know, from LES, the turbulence contributions to the mean Lorentz force and 
compare them with theoretical predictions and DNS results. It is also important to study 
the influence of different effects (e.g., inclination of the initial magnetic field, differential 
rotation, meridional circulations, etc.) on the formation and decay of the inhomogeneous 
magnetic structures. It might be worthwhile to perform LES with larger scale separation, 
and over a broader range of values of the parameters. In this way it should be possible to 
optimize our choice of parameters in order to see where the results of mean-field, DNS, 
and LES can be brought closer together. 
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This study might be important for better understanding of the fundamental mecha- 
nisms governing the formation of sunspots and active regions, and it is expected to have 
an impact in the modelling and predictions of solar magnetic activity. 

We gratefully acknowledge support from Center for Turbulence Research. 
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